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Abstract 

To get a comprehensive view of fungal M35 family (deuterolysin) and M36 family (fungalysin) genes, we conducted 
genome-wide investigations and phylogenetic analyses of genes in these two families from 50 sequenced Ascomycota 
fungi with different life styles. Large variations in the number of M35 family and M36 family genes were found among 
different fungal genomes, indicating that these two gene families have been highly dynamic through fungal evolution. 
Moreover, we found obvious expansions of Meps in two families of Onygenales: Onygenaceae and Arthodermataceae, 
whereas species in family Ajellomycetace did not show expansion of these genes. The strikingly different gene duplication 
and loss patterns in Onygenales may be associated with the different pathogenicity of these species. Interestingly, 
likelihood ratio tests (LRT) of both M35 family and M36 family genes suggested that several branches leading to the 
duplicated genes in dermatophytic and Coccidioides fungi had signatures of positive selection, indicating that the 
duplicated Mep genes have likely diverged functionally to play important roles during the evolution of pathogenicity of 
dermatophytic and Coccidioides fungi. The potentially positively selected residues discovered by our analysis may have 
contributed to the development of new physiological functions of the duplicated Mep genes in dermatophytic fungi and 
Coccidioides species. Our study adds to the current knowledge of the evolution of Meps in fungi and also establishes a 
theoretical foundation for future experimental investigations. 
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Introduction 

Many proteolytic proteases, such as extracellular serine 
proteases, aspartic proteases, and metalloproteases (Meps) are 
known to be important virulence factors related to fungal 
pathogenicity [1,2]. Metalloproteinases, in which zinc is an 
essential metal ion for the catalytic activity, have been identified 
to play important roles in the pathogenicity of pathogenic fungi 
with different life styles, including vertebrate pathogens, entomo- 
pathogens and phytopathogens. For example, a 43.5-kDa Mep 
was identified from the dermatophytic fungus Microsporum canis and 
it could digest keratin azure during infection of cats and guinea 
pigs [3]. A thermolysin-like metalloproteinase produced by the 
entomopathogenic fungus Metarhizium anisopliae has been identified 
to not only cause insect cell and tissue damage [4], but also 
facilitate the development of entomopathogenic fungi within the 
infected insect host [5] . Moreover, a metalloprotease from the rice 
blast fungus Magnaporthe grisea was shown to function as a gene-for- 
gene-specific avirulence factor by directly binding the plant 
resistance gene product and triggering a signaling cascade leading 
to resistance [6] . These studies suggest that Meps have developed 
distinct functional properties to help fungi adapt to different 
ecological niches. Their functional diversities in different fungi 



raise interesting questions regarding the evolution of these 
important genes. 

The Meps identified so far as secreted by pathogenic fungi 
mainly belong to two different families: the deuterolysin (M35) and 
the fungalysin (M36) [2]. These two families of proteases share 
little similarity with each other, with the exception of the HEXXH 
motif, in which two histidine residues function as the first and the 
second zinc ligand, respectively [7]. M35 family genes harbor a 
third zinc ligand known as an Asp that is found in a GTXDXXYG 
motif at the C-terminal, whereas M36 family genes utilize a Glu, 
which is present at the EXXXD motif at the C-terminal as the 
third zinc ligand [7,8]. Previous evidence suggested that members 
of these two families are actual pathogenic factors in human 
infectious diseases. For example, Aspergillus fumigatus, an opportu- 
nistic human pathogen, can release a M36 family metalloprotease 
to help it invade mammalian lung [9] . In addition, at least ten Mep 
genes (designated as Mepl to MeplO) have been found in Coccidioides 
posadasii, which is a fungal respiratory pathogen of humans that 
can cause coccidioidomycosis (Valley fever) in immunocompro- 
mised individuals, and most of them were classified into M35 and 
M36 families. Moreover, the dermatophytic fungi, which cause 
various dermatophytosis common in animals and humans, 
can secrete several M36 family proteases with collagenolytic, 
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elastinolytic and keratinolytic activities to help them survive in 
various niches [10]. Therefore, although the precise roles of these 
proteases from other phytopathogenic or entomopathogenic fungi 
have not yet been identified, we believed that the M35 and M36 
family genes in different pathogenic fungi may have a wide variety 
of pathological actions. Thus, characterizing the evolution of M35 
family and M36 family genes from different fungi can be of great 
importance for understanding the diversification of these Mep 
genes and for inferring their biological significance. Moreover, the 
increasing availability of completely sequenced fungal genomes 
also provides a good opportunity for us to characterize the Meps 
superfamily in different fungal taxa and to infer their evolutionary 
patterns. 

Onygenales, a fungal order associated with a host/substrate shift 
from plants to animals, includes many important pathogenic fungi 
that can cause notorious diseases of vertebrates (e.g. Coccidioides 
immitis, Coccidioides posadasii, Paracoccidioides brasiliensis, Histoplasma 
capsulatum, Trichophyton rubrum and Microsporum canis) and also 
contain a number of non-pathogens (e.g. Uncinocarpus reesii) [11]. 
Previous studies suggested that both the M35 and M36 family 
genes have separately expanded in Onygenales [11,12]. Our 
previous study also provided evidence that positive selection has 
acted on the duplicated M35 family genes in Coccidioides and 
promoted their novel functional adaptations (e.g., recent acquisi- 
tion of virulence traits) of these species [13]. Phylogenetic analyses 
in this paper suggested that M35 and M36 family genes have 
undergone intriguing evolution in Onygenales fungi. To better 
understand the evolutionary force that have likely shaped the 
diversification of these duplicated M35 and M36 family genes in 
Onygenales species and the implications of such changes, we also 
investigated the possible selection pressures responsible for the 
duplicated genes in Onygenales species. 

Data and Methods 

Data sets 

In order to understand the evolutionary history and functional 
divergence of the M35 and M36 family genes in fungi, we 
performed a comprehensive evolutionary analysis of these two 
family genes from the sequenced genomes of 50 Ascomycota fungi. 
These fungi have different life styles and they include saprophytes, 
opportunistic human pathogens, plant pathogens and entomo- 
pathogens. The genome sequences and their corresponding 
annotated proteins database of the 50 Ascomycota fungi are 
available from the Fungal Genome Initiative at the BROAD 
Institute (http://www.broad.mit.edu/annotation/fungi/fgi) or 
Fungal Genome Research (http://fungalgenomes.org/). 

To avoid search bias, we applied two methods to find putative 
homologs of M35 family sequences and M36 family sequences 
from the fungal genome databases employed in this study. In the 
first method, the M35 family genes (Mep2 to Mep8) and M36 
family genes (Mep9 and Mep 10) from C. posadasii were used as 
query sequences to do TBLASTN and BLASTP searches with E- 
value cutoff of 10 -1<) against the genome nucleotide database and 
fungal protein database, respectively. The sequences that met the 
criteria of both searching procedures and with amino acid 
sequence identities greater than 30% and covering more than 
150 amino acids were considered as candidate genes. In the 
second method, the HMM profile Peptidase_M35 (PF02102; 
http://pfam.sanger.ac.uk/ family/PF02 1 02#tabview = tab6) and 
Peptidase_M36 (PF02128; http://pfam.sanger.ac.uk/family/ 
PF02 1 28#tabview = tab6) were downloaded and used as queries 
for homologous protein searches using the program 
HMMSEARCH from the HMMER package (http://hmmer. 



wusti.edu/). For each search result, hits were considered significant 
when they matched the Pfam HMM profile with E values < 10 -5 . 
In total, 105 M35 family genes and 58 M36 family genes were 
identified, including two M35 family genes and four M36 family 
genes from the Basidiomycete fungus Coprinopsis cinerea which were 
used as outgroups for these two families, respectively. The 
taxonomic information and the identified gene number of the 
fungal species are presented in Table 1. 

Phylogenetic analysis 

For each of the two families of genes, MUSCLE v3.5 was used 
to generate protein alignment with default settings [14]. The 
ambiguous areas of alignments were located and removed by using 
the program Gblocks 0.91b [15,16] with default parameters. The 
gap selection criterion "with half was used here. 

Phylogenetic analyses of the alignments were performed using 
PHYML 3.0 [17] for Maximum likelihood (ML) analysis and 
MrBayes 3.1.2 [18] for Bayesian inference. In the ML analysis, the 
best-fitting model of sequence evolution were estimated for both 
M35 family and M36 family genes using program ProtTest [19]. 
The chosen model WAG+I+G and their parameters (I = 0.03, 
G= 1.912 for M35 family genes; 1 = 0.038, G= 1.149 for M36 
family genes) were used in the ML analysis. The reliability of the 
tree topology was evaluated using bootstrap support [20] with 100 
for ML analysis. In addition, Bayesian inference was conducted 
using MrBayes 3.1.2 [18]. The best-fitting model of sequence 
evolution determined by ProtTest [19] was also used as the priors 
of Bayesian inference. Bayesian analysis started with randomly 
generated trees and Metropolis-coupled Markov chain Monte 
Carlo (MCMC) analyses were run for 2 x 10 6 generations. The run 
was stopped when the average standard deviation of split 
frequencies was less than 0.01 in all the cases (MrBayes 3.1.2 
manual). To ensure that these analyses were not trapped in local 
optima, the dataset was run three times independendy. We 
determined the burn-in period by checking for likelihood stability. 
A 50% majority rule consensus of post burn-in trees was 
constructed to summarize posterior probabilities (PPs) for each 
branch. 

Gene duplication and loss analyses 

To investigate the gene duplication and loss scenarios of M35 
and M36 family genes in Onygenales, we used Notung 2.6 [21,22] 
to reconcile the species tree and gene tree. The species tree was 
inferred based on a combined alignment of six genes from each 
Onygenales species as those used in James et al [23], including 18S 
rRNA, 28S rRNA, ITS RNA, translation elongation factor 1-a 
(TEFla), RNA polymerase II largest subunit (RPB1) and RNA 
polymerase II second largest subunit (RPB2) (9601 -bp in total, 
data not shown). Phylogenetic analyses of the alignments were 
performed as described by Li et al [13] and the model TrN+I+G 
of sequence evolution was optimized using Akaike information 
criterion [24,25] as implemented in Modeltest version 3.7 [26]. 
For the gene tree used here, we collapsed those inconsistent nodes 
produced by different tree-building methods, which were also 
poorly supported, into polytomies [27]. 

In addition, we investigated the gene duplication and loss 
scenarios of M35 genes in Eurotiales fungi. The species tree and 
gene tree of Eurotiales fungi were inferred by using the same 
methods mentioned above. 

Selective pressure analyses 

To investigate the possible selective forces behind M35 and 
M36 family gene duplications in Onygenales, we conducted 
likelihood ratio tests (LRT) for these two families respectively. 



PLOS ONE | www.plosone.org 



2 



February 2014 | Volume 9 | Issue 2 | e90225 



Independent Expansion of Meps in Onygenales 



Table 1. Number of both M35 family and M36 family genes in different fungal species. 



Species 


Host 


Taxonomy 


M35 family 


M36 family 


Total 


Coccidioides immitis 


vertebrate pathogen 


Eurotiomycetes:Onygenales:Onygenaceae 


7 


2 


9 


Coccidioides posadasii 


vertebrate pathogen 


Eurotiomycetes:Onygenales:Onygenaceae 


7 


2 


9 


Uncinocarpus reesii 


saprophyte 


Eurotiomycetes:Onygenales:Onygenaceae 


4 


2 


6 


Microsporum canis 


vertebrate pathogen 


Eurotiomycetes:Onygenales:Arthodermataceae 


5 


5 


10 


Microsporum gypseum 


vertebrate pathogen 


Eurotiomycetes:Onygenales:Arthodermataceae 


5 


5 


10 


Trichophyton equinum 


vertebrate pathogen 


Eurotiomycetes:Onygenales:Arthodermataceae 


5 


4 


9 


Trichophyton rubrum 


vertebrate pathogen 


Eurotiomycetes:Onygenales:Arthodermataceae 


5 


4 


9 


Trichophyton tonsurans 


vertebrate pathogen 


Eurotiomycetes:Onygenales:Arthodermataceae 


5 


5 


10 


Paracoccidioides brasiliensis 


vertebrate pathogen 


Eurotiomycetes:Onygenales:Ajellomycetaceae 


1 


- 


1 


Histopiasma capsulatum 


vertebrate pathogen 


Eurotiomycetes:Onygenales:Ajellomycetaceae 


1 




1 


Blastomyces dermatitidis 


vertebrate pathogen 


Eurotiomycetes:Onygenales:Ajellomycetaceae 


2 


- 


2 


Aspergillus flavus 


vertebrate pathogen 


Eurotiomycetes:Eurotiales:Trichocomaceae 


6 


2 


8 


Aspergillus fumigates 


vertebrate pathogen 


Eurotiomycetes:Eurotiales:Trichocomaceae 


3 


1 


4 


Neosartorya fischeri 


vertebrate pathogen 


Eurotiomycetes:Eurotiales:Trichocomaceae 


2 


1 


3 


Aspergillus terreus 


vertebrate pathogen 


Eurotiomycetes:Eurotiales:Trichocomaceae 


2 


1 


3 


Aspergillus clavatus 


saprophyte 


Eurotiomycetes:Eurotiales:Trichocomaceae 


1 


1 


2 


Aspergillus niger 


saprophyte 


Eurotiomycetes:Eurotiales:Trichocomaceae 


- 


1 


1 


Aspergillus nidulans 


saprophyte 


Eurotiomycetes:Eurotiales:Trichocomaceae 


4 


- 


4 


Aspergillus oryzae 


saprophyte 


Eurotiomycetes:Eurotiales:Trichocomaceae 


5 


2 


7 


Penicillium chrysogenum 


saprophyte 


Eurotiomycetes:Eurotiales:Trichocomaceae 


1 


- 


1 


Chaetomium globosum 


endophyte 


Sordariomycetes 


- 


- 


- 


Epichloe festucae 


endophyte 


Sordariomycetes 


- 


- 


- 


Fusarium graminearum 


phytopathogen 


Sordariomycetes 


1 


1 


2 


Fusarium oxysporum lycopersici 


phytopathogen 


Sordariomycetes 


2 


1 


3 


Fusarium solani 


phytopathogen 


Sordariomycetes 


2 


1 


3 


Fusarium verticillioides 


phytopathogen 


Sordariomycetes 


2 


1 


3 


Magnaporthe grisea 


phytopathogen 


Sordariomycetes 


3 


2 


5 


Metarhizium anisopliae 


entomopathogen 


Sordariomycetes 


1 




1 


Neurospora crassa 


saprophyte 


Sordariomycetes 


2 


- 


2 


Podospora anserine 


saprophyte 


Sordariomycetes 


_ 


- 


- 


Verticillium albo-atrum 


phytopathogen 


Sordariomycetes 


4 


2 


6 


Verticillium dahlia 


phytopathogen 


Sordariomycetes 


3 


2 


5 


Trichoderma reesei 


saprophyte 


Sordariomycetes 


1 


1 


2 


Trichoderma atroviride 


saprophyte 


Sordariomycetes 


1 


1 


2 


Trichoderma virens 


saprophyte 


Sordariomycetes 


1 


1 


2 


Phaeosphaeria nodorum 


phytopathogen 


Dothideomycetes 


2 


2 


4 


Pyrenophora tritici-repentis 


phytopathogen 


Dothideomycetes 


1 


1 


2 


Botryotinia fuckeliana 


phytopathogen 


Leotiomycetes 


2 


- 


2 


Botrytis cinerea 


phytopathogen 


Leotiomycetes 


2 


- 


2 


Sclerotinia sclerotiorum 


phytopathogen 


Leotiomycetes 


2 


- 


2 


Candida albicans 


vertebrate pathogen 


Saccharomycotina 


- 


- 


- 


Candida giabrata 


vertebrate pathogen 


Saccharomycotina 




- 


- 


Candida tropicalis 


vertebrate pathogen 


Saccharomycotina 


_ 


_ 


_ 


Debaryomyces hansenii 


saprophyte 


Saccharomycotina 








Pichia stipitis 


saprophyte 


Saccharomycotina 








Saccharomyces cerevisiae 


saprophyte 


Saccharomycotina 








Yarrowia lipolytica 


saprophyte 


Saccharomycotina 








Schizosaccharomyces japonicus 


saprophyte 


Taphrinomycotina 








Schizosaccharomyces pombe 


saprophyte 


Taphrinomycotina 








Coprinopsis cinerea 


saprophyte 


Basidiomycota: Agaricomycotina 


2 


4 


6 



doi:10.1371/journal.pone.0090225.t001 
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From the phylogeny shown in Figure 1, we can see that multiple 
gene duplications of the M35 family genes have occurred in 
Eurotiales, and some of the duplication events occurred after the 
divergence of Eurotiales and Onygenales. To avoid the interfer- 
ence of within-Eurotiales duplications on our analyses, we used 
more distandy related Sordariomycetes fungi as outgroup to 
conduct LRTs analysis of M35 family genes. For M36 family 
genes, we still used Eurotiales fungi as outgroup. Finally, a 861 -bp 
codon-based alignment for 62 M35 family genes (Fig. SI) and a 
1569-bp codon-based alignment for 38 M36 family genes (Fig. S2) 
were used to construct phylogenetic trees using the same methods 
mentioned above except that the best-fitting model of sequence 
evolution used for ML and Bayesian analyses were estimated using 
program Modeltest 3.06 [26]. Because the likelihood analysis 
might be sensitive to tree topology used, we collapsed the nodes 
that showed inconsistent branching patterns from different tree- 
building methods and with poor statistical support into polytomies 
[27]. 

The ratio CO (dN/dS) is the ratio of the number of non- 
synonymous substitutions per non-synonymous site (dN) to the 
number of synonymous substitutions per synonymous site (dS), 
which provides an indication of the change in selective pressures 
[28]. A dN/dS ratio =1, <1, and >1 are indicative of neutral 
evolution, purifying selection, and positive selection on the protein 
involved, respectively [29,30]. The codon substitution models 
implemented in the GODEML program in the PAML 4.4b 
package [31] were used to analyze changes of selective pressure. 
All models were corrected for transition/transversion rate and 
codon usage biases (F3x4). The ambiguous sites were removed in 
PAML analysis (clean data = yes). Different starting CO values were 
also used to avoid the local optima on the likelihood surface [32]. 

Two branch-specific models, i.e., "one-ratio" (MO) and "free- 
ratios", were compared. The MO model assumes the same co ratio 
for all branches while the "free-ratios" model assumes an 
independent co ratio for each branch [33]. We constructed LRT 
to compare the two models. Significant differences between 
models were evaluated by calculating twice the log-likelihood 
difference following a % 2 distribution, with the number of degrees 
of freedom equal to the difference in the numbers of free 
parameters between models. 

Site-specific models, which allow for variable selection patterns 
among amino acid sites, Mia, M2a, M7, and M8, were used to 
test for the presence of sites under positive selection. The M2a and 
M8 models identify positively selected sites. When these 2 positive- 
selection models fitted the data significantly better than the 
corresponding null models (Mia and M7), the presence of sites 
with co> 1 would be suggested. The conservative Empirical Bayes 
approach [34] was then used to calculate the posterior probabil- 
ities of a specific codon site and identify those most likely to be 
under positive selection. 

Considering that positive selection may act in very short 
episodes during the evolution of a protein [35] and affect only a 
few sites along a few lineages in the phylogeny, the likelihood 
models accommodating CO ratios to vary both among lineages of 
interest and among amino acid sites, as an improved version of the 
"branch-site" model [36], were also considered here. We used 
branch-site Model A as a stringency test (test 2) and identified 
amino acid sites under positive selection by an empirical Bayes 
approach along the lineages of interest [36,37]. The log-likelihoods 
for the null and alternative models were used to calculate a 
likelihood ratio test statistic, which was then compared against the 
X2 distribution (with a critical value of 3.84 at a 5% significance 
level) [31]. We used the conservative Bayes Empirical Bayes (BEB) 
approach, which assigns a prior to the model parameters and 



integrates over their uncertainties, to calculate the posterior 
probabilities of a specific codon site and to identify those most 
likely to be under positive selection. In addition, the Bonferroni 
correction [38,39] was also applied for correcting multiple testing 
in the analysis according to the number of tests of significance 
performed. 

Results 

Taxonomic distribution of M35 and M36 gene families 

As shown in Table 1 , large variations of the gene numbers for 
both M35 and M36 family genes were observed from the analyzed 
Ascomycota species, suggesting these two gene families have been 
highly dynamic through fungal evolution. 

For the M35 family, a total of 103 genes were identified from 
the analyzed Ascomycota fungal species. In Family Onygenaceae 
of the Onygenales Order, 7 genes were identified from the human 
pathogenic fungi C. po.sada.sii and C. immitis, and 4 genes were 
predicted from the saprophytic fungus U. reesii. In Family 
Arthodermataceae of Onygenales, five genes each were identified 
from the five dermatophytic fungi (M. cards, Microsporum gypseum, 
Trichophyton equinum, T. rubrum and Trichophyton tonsurans). These 
fungi can cause dermatosis in humans and animals. In Family 
Ajellomycetace of Onygenales, only one M35 family gene each 
was predicted from H. capsulatum and P. brasiliensis, and two genes 
were identified from Blastomyces dermatitidis . This result suggested 
that M35 family genes have expanded in two families of 
Onygenales: Onygenaceae and Arthodermataceae. The genus 
Aspergillus, which includes species with very different life styles with 
some very harmful and others beneficial to humans, showed the 
most variations of gene numbers (range from 0 to 6). For other 
fungal species, most contained 1 to 4 M35 family genes. However, 
no M35 family gene was identified from several fungal groups (e.g. 
Saccharomycotina fungi, Taphrinomycotina fungi, two endophyt- 
ic fungi and some saprophytic fungi) (Table 1). 

For the M36 family, a total of 54 genes were identified from the 
analyzed Ascomycota fungal genome sequences. In Onygenales 
fungi, the dermatophytic fungi (Arthodermataceae) contained the 
most with 4—5 genes per species, while 2 genes each were identified 
from Onygenaceae fungi [C. posadasii, C. immitis and U. reesii) and 
no M36 family gene was identified from Arthodermataceae fungi 
(H. capsulatum, P. brasiliensis and B. dermatitidis), suggesting that M36 
family genes might have expanded in dermatophytic fungi. For 
other fungal species, 1 or 2 M36 family genes were predicated. 
However, several fungal groups (including Saccharomycotina 
fungi, Taphrinomycotina fungi, two endophytic fungi, three 
Leotiomycetes fungi, M. anisopliae, Neurospora crassa and Podospora 
anserine) contain no obvious M36 family gene homologs (Table 1). 
These results suggest that the M36 family genes may have 
undergone many gene duplication and gene loss events during the 
evolution of Ascomycota fungi. 

Phylogenetic analysis of M35 family 

Phylogenetic relationships among the M35 family genes were 
analyzed based on an alignment consisting of 185 amino acids 
from 105 M35 family genes (include 103 genes from Ascomycota 
fungi and 2 outgroup genes from the Basidiomycota fungus C. 
cinerea; Fig. S3) using maximum-likelihood and Bayesian tech- 
niques. As seen in Fig. 1 and Fig. S4, phylogenetic analyses 
consistently grouped the M35 family genes into 12 distinct clades, 
which are abbreviated as M35-a to M35-1. Among these clades, 
M35-a (BS = 99%; PP= 100%) is a dermatophyte-specific lineage, 
containing exclusively the M35 genes from five dermatophytic 
fungi with Microsporum spp. branching basally to Trichophyton spp. It 
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Arthodermataceae 
Onygenaceae 
Eurotiales 
Ajellomycetaceae 
Sordariomycetes 
Leotiomycetes 



Microspomm canis MCYG 04257T0 

'^i — Microspomm gypseum MG7G 02351T0 
U _ Trichophyton rubrum TERG 03599T0 
s** 1 — . Trichophyton tonsurans TE5G_01292T0 
n Trichophyton equinum TEQG 08048T0 
. Microspomm cams MCYGJI991T0 

- Microspomm gypseum MGYG_04094T0 
. Trichophyton rubrum TERG 04416T0 
. Trichophyton tonsurans TESG 06846T0 
I Trichophyton equinum TEQG ~06169T9 

. Microspomm canis MCYG J152Q1T0 

. Microspomm gypseum MGYG 0624IT8 

. Trichophyton rubrum TERG 06377T0 

' 1 , Trichophyton tonsurans TESG_03882T0 

I. Trichophyton equinum TEQG_0467IT0 
Microspomm canis MCYG 07IS5T0 
Microspomm gypseum MGYG 03465T0 
Trichophyton rubrum TERGJ1528T0 

— Trichophyton tonsurans TESG 05I02T0 
Trichophyton equinum TEQG059Z2T0 

Microspomm gypseum MGYG 608I3T0 
Microsporia canis MCYG~ 00239T9 
Trichophyton rubrum TERG 00673T0 
Trichophyton tonsurans TESG 02193T0 
Tri chophyton eauinum TEQG 02633T0 




M35-a 



,. Coccidhides posadasU CPA T 05050 
H Coccidioides immitis CIMG S86U.ll 
. Coccidhides posadasU CP.fT 07671 
i Coccidioides immitis CIMG TOlOl.tl 
Coccidhides posadasU CPAT 04742 
Coccidhides immitis CIMG 07349.11 



M3S-b 



Uncinocarpiis reesii LJRET 01255 

Uncmocarpus reesiiVRET 02006 
Unciiwcarpits reesii V RET 03761 
Coccidhides posadasU CPAT 04075 
Coccidioides immitis CIMG 06682.11 
Coccidhides posadasU CPAT 08585 
Coccidhides immitis CIMG 05736.11 
Coccidhides posadasU CPAT 08667 
Coccidioides immitis CIMG OlOlO.tl 



M35-C 



Seosartona jischeri ffl'H 1026J0.ll 

Aspergillus fumigates A FUAJG1J 750.11 
, Neosartorya fisclteri NFIA 031120.11 
11 Aspergillus fumigates AFVA 4G02700.ll 
Aspergillus clavalus A CIA 052 720.il 
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Fusarium verticillioides FVET 01806 
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Figure 1. ML tree based on amino acid sequences of 105 M35 family genes. The tree was performed using PHYML 3.0[17]. The best-fitting 
model WAG+I+G and their parameters (I = 0.03, G = 1 .91 2) which were estimated by program ProtTest [50] were used in the ML analysis. The reliability 
of the tree topology was evaluated using bootstrap support [20] with 100. 
doi:1 0.1 371 /journal.pone.0090225.g001 



seemed that the multiplication of these genes occurred before the 
divergence of these five species. This is because each subclade 
within this clade contained one gene from each species, suggesting 
that these genes had a conserved and essential function in 
dermatophytes. The genes from Onygenaceae fungi fell into three 
clades (M35-b, M35-c and M35-d). Clades M35-c (BS and 
PP= 100%) and M35-d (BS and PP= 100%) were comprised of 
sequences from C. posadasii, C. immitis and U. reesii, while M35-b 
(BS = 88%; PP = 98%) was a Coccidioides-specific lineage, consisted 
of three duplicated genes from each of the two Coccidioides species. 
The genes from Ajellomycetaceae fell into clades M35-e (BS and 
PP = 100%) and M35-J (BS and PP= 100%). Clade M35-J, which 
contained genes from B. dermatitidis and P. brasiliensis, was closely 
related to Onygenaceae lineage M35-d. In contrast, clade M35-e, 
which included one gene from H. capsulatum and another from B. 
dermatitidis, showed close relationships with Eurotiales lineage 
M35-g (BS = 92%; PP = 1 00%). 

From our phylogenetic analysis, we found that except for 
Arthodermataceae and Onygenaceae, the M35 family genes from 
other Ascomycota fungi within the same categories were 
consistently clustered into two clades respectively. For example, 
the genes from Ajellomycetaceae fell into M35-e and M35-J clades 
and those from Eurotiales were clustered into M35-g and M35-h 
clades (BS = 89%; PP = 86%). Clades M35-i (BS and PP= 100%) 
and M35-1 (BS = 67%; PP = 70%) were comprised of genes from 
Sordariomycetes while clades M35-j (BS and PP = 100%) and M3- 
k (BS and PP= 100%) mainly consisted of those from Leotiomy- 
cetes fungi. Noticeably, one gene from Meurospora crassa showed a 
close relationship with clade M35-j. Interestingly, although the 
genes from the same taxonomic groups of fungi were consistendy 
clustered into two clades, gene duplication and loss events might 
have also occurred during their evolution, leading to the 
diversification of M35 family genes in different fungi. For instance, 
though clades M35-i and M35-1 both contained genes from 
Sordariomycetes, clade M35-1 contained more genes than clade 
M35-L 

Phylogenetic analysis of M36 family 

Results from our phylogenetic analyses based on an alignment 
consisting of 538 amino acids from 58 M36 family sequences 
(include 54 genes from Ascomycota fungi and 4 outgroup genes 
from the Basidiomycota fungus C. cinerea, Fig. S5) were largely 
consistent with the taxonomical relationships among the major 
fungal groups. The results suggested that the Dothideomycetes 
showed close relationships with Sordariomycetes and that the 
Eurotiales species were closely related to Onygenales fungi. 
Phylogenetic analyses based on the M36 family genes consistendy 
showed eight strongly-supported monophyletic clades which were 
designated as M36-1 to M36-8 in the trees (see Figure 2 and Fig. 
S6). The M36 family genes from Arthodermataceae fungi (M. 
canis, M. gypseum, T. equinum, T. rubrum and T. tonsurans) and 
Onygenaceae fungi fell into two clades, respectively. Among these 
clades, clades M36-1 (BS and PP = 100%) and M36-3 (BS = 99%; 
PP= 100%) each contained the M36 genes from dermatophytic 
fungi. These genes seemed to have duplicated before the 
divergence of these five dermatophytic species. The sequences 
from Onygenaceae fungi (C. posadasii, C. immitis and U. reesii) 
grouped into clades M36-2 (BS = 94%; PP= 100%) and M36-4 
(BS and PP = 100%), and their duplication likely occurred after the 



divergence of Onygenaceae and Arthodermataceae families. The 
genes from Eurotiales and Dothideomycetes fell into clades M36-5 
(BS = 77%; PP = 94%) and M36-6 (BS = 68%; PP = 99%), respec- 
tively. The genes from Sordariomycetes were mainly grouped into 
clades M36-7 and M36-8 (BS and PP= 100%). However, clade 
M36-8 was comprised of the duplicated genes from Verticillium albo- 
atrum and Verticillium dahlia, showing distant divergence with other 
clades. 

Gene duplication and loss analysis in Onygenales and 
Eurotiales 

Based on the above results, we found that M35 and M36 family 
genes have undergone expansion in Onygenales fungi. To infer 
the gene duplication and loss scenarios of M35 and M36 family 
genes in Onygenales, We used Notung 2.6 [21,22] to reconcile the 
species tree (Fig. S7) and the gene tree of each family. Using 
various tree-building methods, our inferred species tree consis- 
tently showed three strongly supported clades (Fig. S3), consistent 
with previous reports, which suggested that the fungal family 
Onygenaceae was more closely related to Arthodermataceae than 
to Ajellomycetace [13,40-42]. The inferences of gene duplication 
and loss events in Onygenales species are shown in Figure 3. 

For the M35 family genes (Figure 3a), there were at least 10 
gene duplication and 5 gene loss events in the Onygenales species. 
An initial duplication event occurred before the divergence of 
Onygenales fungi. Subsequendy, four and five gene duplications 
occurred independently in Arthodermataceae and Onygenaceae, 
respectively. In Onygenaceae, two gene duplications occurred 
before the divergence of Coccidioides and U. reesii, while three 
occurred before the divergence of C. posadasii and C. immitis 
(Figure 3a). For comparison, one M35 gene might have been lost 
in each of U. reesii, P. brasiliensis and H. capsulatum. Moreover, the 
lineages of Arthodermataceae might have also suffered two gene 
loss events after the duplication but before the divergence of the 
five dermatophytic fungi (Figure 3 a). 

For the M36 family genes (Figure 3b), there were four gene 
duplication and four gene loss events in Onygenales species. 
Before the divergence of Onygenales, there was one initial 
duplication event. Subsequendy, three duplication events occurred 
in Arthodermataceae. There was one gene loss event in each of the 
two lineages of Onygenaceae and Arthodermataceae respectively 
and one gene was also lost in each of T. equinum and T. rubrum. 

In our study, we also found large variations of M35 family gene 
number in Eurotiales, ranging from 0 to 6 (Table 1). To infer the 
gene duplication and loss scenarios of M35 genes in this group, we 
also used Notung 2.6 [21,22] to reconcile the species tree (Fig. S8) 
and the gene tree of this family. The inferences of gene duplication 
and loss events in Eurotiales species are shown in Figure 4. We 
found at least 7 gene duplication and 14 gene loss events in the 
Eurotiales species. While most of the duplication events occurred 
before the divergence of Aspergillus spp, the gene loss events in 
Eurotiales fungi were more dynamic. There were three gene loss 
events occurred in Aspergillus nidulans and Aspergillus terreus each, and 
one gene might have been lost in Aspergillus oryzae, Aspergillus clavatus 
and Penicillium chrysogenum each. Moreover, four gene loss events 
occurred in the lineage including three species {A. Jumigatus, 
Neosartorya Jischeri and A. clavatus) and one gene might have been lost 
before the divergence of the four species including A. nidulans, A. 
terreus, A. oryza and Aspergillus jlavus. 
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Arthodermataceae 

Onygenaceae 

Eurotiales 

Dothideomycetes 

Sordariomycetes 




75 



72 



SO 



Trichophyton equinum TEQG 071471T0 
Trichophyton tonsurans TESGJ7429T0 
Trichophyton rubrum TERGJ2213T0 
Mkrosporum gypseum MGYG 05327T0 
Microsporum canis MCYG 07415T0 

I Trichophyton equinum TEQGJ0173T0 

I ' Trichophyton tonsurans TESGJ3496T0 

' Microsporum gypseum MGYG 05070T0 

Microsporum canis MCYG_0768ST0 



M36-1 



Uncinocarpus reesii URET 06232 
Coccidioides posadasii CPAT 03439 
Coccidioides immitis CIMG 06073.tl 



M36-2 



99 



ma 



93 



L 

,-P- Tri, 
0T\— Mil 

I Mir 



Trichophyton equinum TEQG 04323T0 
Trichophyton tonsurans TESG 0040IT0 
Trichophyton rubrum TERG 04324T0 
Microsporum gypseum MGYG 03984T0 

Microsporum canis MCYG 02088T0 

Trichophyton equinum TEQG 05259T0 
Trichophyton tonsurans TESG 06485T0 
Trichophyton rubrum TERG 04809TO 
Microsporum gypseum MGYG 02992T0 
Microsporum canis MCYGJ2168T0 
Trichophyton rubrum TERG_03248T0 
Trichophyton tonsurans TESGJI277ST0 
Microsporum gypseum MGYG 08200T0 
Microsporum canis MCYG 04620T0 



M36-3 



100 



Coccidioides immitis CIMG_1019I.tl 
Coccidioides posadasii CPA T 07 7 63 
— Uncinocarpus reesii URET 03538 



M36-4 



11 



_ r A s 



c 



Aspergillus terreus ATEG_07S44.il 

Aspergillus niger An01g02070.tl 

Aspergillus oryzae AO090I3S000160.il 
Aspergillus flavus AFL2T 08800 
Aspergillus clavatus ACLA_060300.ll 
Neosartorya fisclieri NF1A 099860.H 
Aspergillus fumigates AFUA8G07080.il 

r Aspergillus oryzae AO0900U000036.il 

I Aspergillus flavus A FL2T 04842 



M36-5 



68 



Phaeosphaeria nodorum Pnod SN15-.SNOG 00493.U 

Pyrenophora tritici-repentis PTRG 10639.U 

Phaeosphaeria nodorum Pnod_SN15:SNOGJ0695.tl 



M36-6 



100 



Magnaportlte grisea MGG 13250.U 
Magnaporthe grisea MGG 13252. 11 

Trichoderma atroviride TRIATDRAFT 299629 



100 



Trichoderma virens TRIVIDRAFT 85599 

— Trichoderma reesei tree_ 70800 

— Fusarium oxysporum lycopersici FOXT16612 

Fusarium verticillioides FVET 13650 

Fusarium graminearum FGST 04568 

Fusarium solani fso/ 8789 



M36-7 



Verticillium dahlia VDAG 03418T0 



Verticiliium albo-atrum VDBG 06264T0 



100 



CI 



Verticillium dahlia VDAG 01216T0 



100 



Verticillium albo-atrum VDBG 03987T0 

• Coprinopsis cinerea CC1G_11771T0 
Coprinopsis cinerea CC1G 12248T0 
— Coprinopsis cinerea CC1G 09223T0 
Coprinopsis cinerea CC1G_14256T0 



M36-8 



Outgroup 



0.2 

Figure 2. ML tree based on amino acid sequences of 58 M36 family genes. The tree was performed using PHYML 3.0[17]. The best-fitting 
model WAG+I+G and their parameters (1 = 0.038, G = 1.149) which were estimated by program ProtTest [50] were used in the ML analysis. The 
reliability of the tree topology was evaluated using bootstrap support [20] with 100. 
doi:1 0.1 371 /journal.pone.0090225.g002 



Selective pressure analyses in Onygenales 

Table 2 shows the evidence of positive selection analyses of M35 
family genes in Onygenales. In the branch-specific model analyses 



[43], the free-ratio model, Ml, revealed a significandy better fit to 
the data than did the one-ratio model, M0 (p<0.001, Table 3), 
suggesting that these M35 family genes have been the subjects of 
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Microsporum canis MCYG 042S7T0 
Microsporum gypseum MGYG 0235 ITU 
Trichophyton rubrum TERG H3599T0 
Trichophyton tonsurans TESG_01292T0 
Trichophyton equinum TEQG 0S048T0 
Microsporum canis MCYG 05291 TO 
Microsporum gypseum MGYG 06241T0 
Trichophyton rubrum TERG 0S377T0 
Trichophyton tonsurans TESGJI3882TII 
Trichophyton equinum TEQG 0467 1 TO 
Microsporum canis MCYG 0T155T0 
Microsporum gypseum MGYG 03465T0 
Trichophyton rubrum TERG 01528T0 
Trichophyton tonsurans TESG 05I02T0 
Trichophyton equinum TEQG 05962T0 
Microsporum canis MCYG 0T991T0 
Microsporum gypseum MGYG 04094T0 
Trichophyton rubrum TERG U4416T0 
Trichophyton tonsurans TESG_06846T0 
Trichophyton equinum TEQG 06I69T0 
Microsporum canis MCYG 00239T0 
Microsporum gypseum MGYG 00813T0 
Trichophyton rubrum TERG IT0673T0 
Trichophyton tonsurans TESG 02I93T0 
Trichophyton equinum TEQG 02633T0 
Coccidioides posadasii CPAT 05050 
Coccidioides immitis CIMG U8613.tl 
Coccidioides posadasii CPAT 07671 
Coccidioides immitis CIMG TOlOl.tl 
Coccidioides posadasii CPAT 04742 
Coccidioides immitis CIMGJT7349.fi 
Vncinocarpus reesii LOST* 
Coccidioides immitis CIMG 03010AI 
Coccidioides posadasii CPAT 08667 
Coccidioides immitis CIMG 05736.fl 
Coccidioides posadasii CPAT 08585 
Vncinocarpus reesii URET 01255 
Coccidioides immitis CIMG 06682.fl 
Coccidioides posadasii CPAT 04075 
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Blastomyces dermutitidis BDCG 00922 
Histoplasma capsulatum LOST*~ 
Paracoccidioides brasiliensis PADG_00776T0 
Coccidioides immitis CIMG 00508 
Coccidioides posadasii CPAT_02396 
Vncinocarpus reesii URET_04198 
Arthrlodermatuceae LOST* 
Blastomyces dermutitidis BDCG 03454 
Histoplasma capsulatum HCAG~_05788.tl 
Paracoccidioides brasiliensis LOST* 



Microsporum canis MCYG_07688T0 
Microsporum gypseum MGYG_05070T0 
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Trichophyton tonsurans TESG 03496T0 
Trichophyton equinum TEQG00I73T0 
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Trichophyton rubrum TERGJ2213T0 
Trichophyton tonsurans TESG_07429T0 
Trichophyton equinum TEQG 07147IT0 
Vncinocarpus reesii VRET 06232 
Coccidioides posadasii CPA T 03439 
Coccidioides immitis CIMG_06073.tl 
Trichophyton equinum TEQG 04323T0 
Trichophyton tonsurans TESG 0040IT0 
Trichophyton rubrum TERG 04324T0 
Microsporum gypseum MGYG 03984T0 
Microsporum canis MCYGJ2088T0 
Trichophyton tonsurans TESG06485T0 
Trichophyton equinum TEQG 05259T0 
Trichophyton rubrum TERGJ4809T0 
Microsporum gypseum MGYG_02992T0 
Microsporum canis MCYG 02168T0 
Trichophyton equinum LOST* 
Trichophyton tonsurans TESG_02775T0 
Trichophyton rubrum TERGJ3248T0 
Microsporum gypseum MGYG 08200T0 
Microsporum canis MCYG 04620T0 
Onygenaceae LOST* 
Coccidioides immitis CIMG_10I9I.tl 
Coccidioides posadasii CPAT 07763 
Vncinocarpus reesii VRET 03538 
Arthrlodermataceae LOST* 
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Figure 3. Duplication and loss events of M35 family genes and M36 family genes in Onygenales fungi. The reconciliation between 
species tree and gene tree of Onygenales fungi along with the confirmation of the gene loss/duplication scenario were determined by using Notung 
2.6 [21]. The species tree of Onygenales fungi is shown as Fig. S7. Putative duplication events are indicated with solid cycles, while loss events are 
indicated with thick branches, a, duplication and loss events of M35 family genes, b, duplication and loss events of M36 family genes. 
doi:1 0.1 371 /journal.pone.0090225.g003 

different selective pressures. In the site-specific model analyses, significant in four of the branches (branch a, e, f and j; 
although the LRT of M2a/Mla did not achieve statistical p<0. 002632) (Table 2, Figure 5). Remarkably, several positively 
significance (P= 1.000), M8, another positive-selection model selected residues (A150R/Q, K166A/R/G and M203R/H for 
provided a significantly better fit to the data than did the neutral branch a, V177Q and Y255Q, for branch e, E57T, G79K/Q, 
model (M7) (P<0.001), suggesting the possibility of positive E85S/T/L and A99K for branch /; Y63V/I and A84R for branch 
selection acting on the M35 family genes in Onygenales examined j) were also identified for these branches with high posterior 
here. When we performed the branch-site model tests for those probabilities (Table 2 and Figure 5). Among these eleven positive 
branches resulted from gene duplications of M35 family in selected sites, ten of them (A150R/Q, K166A/R/G, M203R/H, 
Onygenales (19 branches in total, a-t as indicated in Figure 5), we V177Q, E57T, G79K7Q, E85S/T/L, A99K, Y63V/I and A84R) 
found that there were five branches (branches a, d, e, f and j) were involved in amino acids substitutions with different physico- 
showing signs of positive selection (Figure 5). After Bonferroni chemical properties (Table. SI), 
correction for multiple testing, we found that LRT tests were still 

Neosartorya fischeri NFIA_102630.il 
Aspergillus fumigates AFUA_4G13750.tl 
Neosartorya fischeri NFIA 031120.H 
Aspergillus fumigates AFUA_4G02700.il 
Aspergillus clavatus ACLA_052720.tl 
Aspergillus oryzae AO090001000l3S.il 
Aspergillus flavus AFL2T 07373 
Aspergillus terreus ATFG_04941.ll 
Aspergillus nidulans AN7962.2.I1 
Aspergillus nidulans AN3393.2.I1 
Aspergillus oryzae AO090010000493.il 
Aspergillus flavus AFL2T11655 
Aspergillus terreus LOST * 
1 LOST * 

Aspergillus nidulans LOST * 
Penicillium chrysogenum Pc06g00080.tl 
Aspergillus oryzae AO090701000313.il 
Aspergillus flavus AFL2T_05945 
Aspergillus terreus ATEG_08109.il 
1 LOST * 

Aspergillus nidulans AN3959.2.U 
Aspergillus oryzae AO09000500l350.il 
Aspergillus flavus AFL2TJ1268 
Aspergillus flavus AFL2TJ8988 
Aspergillus oryzae LOST * 
Aspergillus terreus LOST * 
1 LOST * 

Aspergillus nidulans LOST * 
Aspergillus oryzae AO090038000111.il 
Aspergillus flavus AFL2T 08935 
Aspergillus terreus LOST * 

1 LOST * 

Aspergillus nidulans LOST * 
Aspergillus nidulans AN83I2.2.H 
Aspergillus fumigates AFUA_2G18070.il 
Neosartorya fischeri LOST * 
Aspergillus clavatus LOST * 

2 LOST * 

Penicillium chrysogenum LOST * 

Figure 4. Duplication and loss events of M35 family genes in Eurotiales fungi. The reconciliation between species tree and gene tree of 
Eurotiales fungi along with the confirmation of the gene loss/duplication scenario were determined by using Notung 2.6 [21]. The species tree of 
Eurotiale fungi is shown as Fig. S8. Putative duplication events are indicated with solid cycles, while loss events are indicated with thick branches. 1, 
the lineage includes species of A. fumigatus, N. fischeri and A. clavatus. 2, the lineage includes species of A. nidulans, A. terreus, A. oryza and A. flavus. 
doi:1 0.1 371 /journal. pone.0090225.g004 
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Figure 5. Phylogenetic tree of 62 M35 family genes used for codon-based maximum likelihood analysis in PAML. Phylogenetic trees 
were collapsed with inconsistent nodes from different tree-building methods and poor statistical supports into polytomy. Branches a-s indicated 
putative duplication events in Onygenales fungi. The branches with significant evidence of positive selection are indicated as a thick branch. The 
putative positively selected residues along these branches were shaded in grey. 
doi:1 0.1 371 /journal.pone.0090225.g005 



Table 3 shows the results of positive selection analyses of M36 
family genes in Onygenales. In the branch-specific model analyses 
[43], the free-ratio model, Ml, revealed a significantly better fit to 
the data than did the one-ratio model, MO (p<0.001, Table 3), 
suggesting that these M36 family genes have been the subjects of 
different selective pressures. Similar as the results in M35 family, 
only the LRT of M8/M7 achieved statistical significance 
(P<0.001) in the site-specific model analyses, also suggesting the 
possibility of positive selection acting on the M36 family genes in 
Onygenales examined here. Intriguingly, the LRT tests based on 
the branch-site models for those branches resulted from gene 
duplications of M36 family in Onygenales (12 branches in total, a-l 
as indicated in Figure 6) suggested that there was significant 
evidence along several lineages leading to the duplicated genes in 
dermatophytic fungi (branch a, b and g). After performing 
Bonferroni correction for multiple testing, the LRT tests still 
showed significant results in these branches (p<0.005). In 
addition, several residues (G183P, Y267V and K401H for branch 



a, T40E for branch b, A139H, I221A/P, Y232W, L300C and 
W510M for branch g) were predicted as positively selected sites 
along these branches (Table 3 and Figure 6). Among these nine 
positive selected sites, five of them (G183P, Y267V, T40E, 
A139H, Y232W) were involved in amino acid substitutions with 
different physico-chemical properties (Table. SI) 

Discussions 

To get a comprehensive view of fungal M35 family and M36 
family genes, we conducted genome-wide investigations and 
phylogenetic analyses of these two family genes from 50 sequenced 
Ascomycota fungal genomes with different life styles. Our results 
provided valuable insights for understanding the evolution of these 
two gene families in various fungi. In Sordariomycetes, the 
saprophytic fungi had 0 to 2 Mep genes, whereas the phytopath- 
ogenic fungi contained 2 to 6 genes within each species. 
Phylogenetic analyses suggested that the M35 family genes of 
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Table 2. CODEML analyses of selective pattern for M35 family genes. 





Models 




\nf 


Parameter Estimates 


2AL b 


Positively Selected Sites' 


Branch-specific models 


MO 


-16572.48016 


til = 0.25237 


455.838*** 






Mia 


-16344.56115 








Site-specific models 


Mia 


— 16214.65883 


0)0 = 0. 1 7007,0)1 = 1 ,p0 = 0.67270, 
p1 =0.32730 




Not allowed 




M2a 


-16214.65883 


0)0 = 0. 1 7007,0)1 = 1 ,0)2 = 1 ,p0 = 0.67270, 
pi =0.23742,p2 = 0.08988 




None 




M7 


- 1 5986.8697 


p = 0.63889,q = 1 .50238 


425.426*** 


Not allowed 




M8 


— 16199.58278 


p = 0.251 1 7,q = 1 .45893,p0 = 1 .00000, 
p1 =0.00000,0) = 2.85555 




55(0.999), 76(0.983) 


Branch-site Branch a 
models 


Null 


-29907.40643 


0)0 = 0. 1 9365,0)1 = 1 ,0)2 = 1 ,p0 = 0.60899, 
p2a = 0.03351 ,p2b = 0.01 865 


12.5344*** 


150 (0.993), 

1 66(0.991 ),203(0.995) 






Alterative 


-29901.13922 

0)0 = 0. 1 9364,0)1 = 1 ,0)2 = 36.77462,p0 = 0.60444, 
p2a = 0.03748,p2b = 0.02091 






Branch e 


Null 


-29906.27165 


0)0 = 0. 1 9306,0)1 = 1 ,0)2 = 1 ,p0 = 0.59293, 
p2a = 0.04993,p2b = 0.02774 


30.5388*** 


177(0.998),255(0.997) 






Alterative 


-29891.00224 

0)0 = 0.19391,0)1 = 1,0)2 = 44.75155, 

pO = 0.55743,p2a = 0.08606,p2b = 0.04768 






Branch f 


Null 


-29901.87988 


0)0 = 0. 1 9307,0)1 = 1 ,0)2 = 1 ,p0 = 0.47997, 
p2a = 0.17002,p2b= 0.09155 


37.969704*** 


57(0.995),79(0.997),85 
(0.997),99(0.991) 






Alterative 


-29882.89503 

0)0 = 0. 1 9340,0)1 = 1 ,0)2 = 98.1 51 94,p0 = 0.553 1 1 , 
p2a = 0.09909,p2b = 0.05284 






Branch j 


Null 


-29907.74419 


0)0 = 0. 1 9363,0)1 = 1 ,0)2 = 1 ,p0 = 0.55036, 
p2a = 0.0943 1,p2b= 0.05198 


9.352552*** 


63(0.983), 84(0.980) 






Alterative 


-29903.06791 

0)0 = 0. 1 9389,0)1 = 1 ,0)2 = 21 .5 1 283,p0 = 0.61 895, 
p2a = 0.02686,p2b = 0.01 473 







a \nL is the log-likelihood scores. 

b LRT to detect adaptive evolution. *** P-C0.001 

Posterior probabilities value of each codon site were showed in parentheses. 
doi:1 0.1 371 /journal.pone.0090225.t002 



Sordariomycetes fungi fell into two clades: M35-i and M35-1. 
Clade M35-1 contained the genes from Sordariomycetes that 
included saprophytic, phytopathogenic and entomopathogenic 
fungi, while clade M35-I comprised six duplicated genes from six 
phytopathogenic fungi each (Figure 1). In addition, the M36 
family genes in Sordariomycetes fungi also grouped into two 
clades: M36-7 and M36-8, with M36-8 containing the duplicated 
genes from the two phytopathogenic fungi V. albo-atrum and V. 
dahlia (Figure 2). Similarly, M35 family genes from the three 
Leotomycetes fungi [Botryotinia fuckeliana, Botrytis cinerea and 
Sclerotinia selerotiorum) were also clustered into two clades: M35-j 
and M35-k, indicating that the genes have duplicated before the 
divergence of these three phytopathogenic fungi. Previous studies 
suggested that protease gene family expansion appeared to be an 
important evolutionary step among the phytopathogenic fungi and 
likely served important roles in plant infection [44] . We speculate 
that the duplication of Mep genes occurred in the phytopathogenic 
fungi belonging to Sordariomycetes and Leotomycetes may be 
associated with their pathogenicity. However, we did not identify 
any Mep gene from the two endophytic fungi, Chaetomium globosum 
and Epichloe festucae, and no M36 family gene was identified from 
the three Leotomycetes fungi, suggesting that these genes may be 
lost during the evolution of these fungi. 



Based on the homologous search results, large variations of gene 
number, ranging from 1 to 8 were also observed in Eurotiales 
(Table 1). Phylogenetic analyses and gene duplication/loss 
analyses suggested that the Mep genes, especially the M35 family 
genes have been highly dynamic through the evolution of 
Eurotiales fungi (Figure 1, 2 and 4). In our study, A. oryzae, which 
is a saprophytic fungus and plays a key role in the fermentation 
process of several traditional Japanese beverages and sauces [45] , 
contains 7 Meps in total (5 in M35 family and 2 in M36 family), 
whereas A. fumigatus, an important opportunistic pathogen of 
vertebrates with compromised immune systems [46], has only 4 
Mep genes in total (3 in M35 family and 1 in M36 family). 
However, although A. fumigatus contains fewer Mep genes, it has 
been reported that this fungus can utilize the M36 family 
metalloprotease to help it invade mammalian lung and degrade 
its host's proteinaceous structural barriers [9]. Moreover, other 
metalloproteases, such as Aspf2 from A. fumigatus and Aspndl from 
A. nidulans which are known cell-surface antigens during fungal 
infections and can bind specific ligands within mammalian hosts, 
were also previously identified from Aspergillus spp. [47]. Therefore, 
we believed that the opportunistic human pathogens in Eurotiales 
might not only use the M35 or M36 family genes as virulence 
factors, but also choose other families of metalloproteases as 
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Table 3. CODEML analyses of selective pattern for M36 family genes. 



Models 




ln£ a 


Parameter Estimates 


2AL" 


Positively Selected Sites' 


Branch-specific models 


MO 


-20322.45718 


co = 0.1 4200 


348.64117*** 






M1 a 


— 20148.1 366 








Site-specific models 


Mia 


- 1 9897.27048 


0)0 = 0.09291,0)1 =1,p0= 0.80509,p1 =0.19491 




Not allowed 




M2a 


- 1 9897.27048 


0)0 = 0.09291,0)1 =1,O)2 = 1,p0 = 0.80509,pl =0.13245, 
p2 = 0.06245 




None 




M7 


-19953.94322 


p = 0.38777,q = 1 .85670 


472.167858** 


f Not allowed 




M8 


-20190.02715 


p = 1 .04901 ,q = 1 .52775,p0 = 1 .00000, 
pi =0.00000,0 = 2.59980 




158 (0.999), 172(0.999), 
198(0.992), 208 (0.992), 

tci m nnoi 'Don in nno 

251(0.998), 380 (0.996), 
391 (0.993) 


Branch-site Branch a 
models 


Null 


-19892.13564 


0)0 = 0.091 67,0)1 = 1 ,0)2 = 1 ,p0 = 0.72573, 
p2a = 0.08330,p2b = 0.01 966 


9.52853*** 


183 (0.996), 
267(0.996),401(0.993) 






Alterative 


-19887.37138 

0)0 = 0.09239,0)1 = 1 ,0)2 = 30.85099,p0 = 0.791 1 4, 
p2a = 0.01 902,p2b = 0.00446 






Branch b 


Null 


- 1 9890.9925 


0)0 = 0.091 54,0)1 = 1 ,0)2 = 1 ,p0 = 0.73420,p2a = 0.07504, 
p2b= 0.01769 


13.15426*** 


40(0.996) 






AltGrstivG 


— 19884.41539 

0)0 = 0.0921 2,0)1 = 1 ,0)2 = 1 5.1 0889,p0 = 0.781 20, 
p2a = 0.03033,p2b = 0.00704 






Branch g 


Null 


-19884.48335 


0)0 = 0.09220,0)1 = 1 ,0)2 = 1 ,p0 = 0.68433, 
p2a = 0.13556,p2b= 0.02978 


21.171318*** 


1 39(0.994), 






Alterative 


221(0.995),232(0.996),300(0.990)510(0.998) 
-19873.89769 

0)0 = 0.09356,0)1 = 1 ,0)2 = 1 1 .57996,p0 = 0.76788, 
p2a = 0.05883,p2b = 0.01 233 







a lni. is the log-likelihood scores. 

b LRT to detect adaptive evolution. *** P-C0.001 

Posterior probabilities value of each codon site were showed in parentheses. 
doi:1 0.1 371 /journal.pone.0090225.t003 



pathogenic factors to carry out the infection of immunocompro- 
mised individuals. 

Surprisingly, neither M35 gene nor M36 gene was found in 
Saccharomycotina and Taphrinomycotina fungi, suggesting that 
both the two family genes have been lost during the evolution of 
the yeast-like fungi. However, although no M35 or M36 family 
proteins have been identified in yeast-like fungi, other types of 
metalloproteases that serve important roles in the development 
and pathogenicity of yeast-like fungi have been identified before 
[48]. For example, a Zpslp metalloproteases which involved in 
zinc-responsive transcriptional regulation was observed in Saccha- 
romyces cerevisiae [49]. Moreover, a Zpslp-like metalloproteases 
named Pral which function as virulence factors in mediating 
critical host-parasite interactions were identified from the oppor- 
tunistic human pathogen Candida albicans that can cause diseases of 
vertebrates [48] . Sequence alignments suggested that some critical 
residues such as the cysteine residues involved in disulfide bond 
formation and a zinc-binding domain were conserved in both 
Zpslp and M35 family metalloproteases [48]. These results 
suggest that the yeast-like fungi may utilize other metalloproteases 
but not M35 or M36 family proteins to carry out their 
development and pathogenicity. 

The most striking result from our study is the significant 
varieties of Meps in Onygenales fungi. Phylogenetic reconstruction 
and gene duplication/loss analyses suggest that the first duplica- 
tion event of both M35 and M36 family genes occurred in the 
common ancestor of the Onygenales species (Figure 1, 2 and 3). 



Subsequently, additional gene duplication events mainly occurred 
after the divergence of Arthrlodermataceae (dermatophyte) and 
Onygenaceae species (Figure 1, 2 and 3), leading to a large 
number of Mep genes shared between Arthrlodermataceae 
(dermatophyte) and Onygenaceae species. However, the Ajello- 
mycetace species (B. dermatitidis, H. capsulatum and P. brasiliensis) did 
not show evidence of expansion of these genes (only containing 1 
or 2 M35 family genes and no M36 family gene was identified 
from them). Although the species in Onygenales showed a host/ 
substrate shift from plants to animals during evolution, their 
pathogenicity factors on animals might be different. Dermato- 
phytic fungi cause the majority of superficial dermatophytosis in 
humans and animals [12], while Coccidioides spp. can cause life- 
threatening respiratory disease known as coccidioidomycosis 
(Valley fever) in immunocompromised individuals [50,51]. How- 
ever, Ajellomycetace species cause geographically widespread 
systemic mycosis of humans and other mammals [52]. For 
example, Paracoccidioides can cause paracoccidioidomycosis [40], 
B. dermatitidis is the causative agent of blastomycosis. However, H. 
capsulatum is the most common cause of histoplasmosis in the world 
[53]. Therefore, the distinct expansion of Meps in Arthrloderma- 
taceae and Onygenaceae species suggest that the Meps may have 
played different roles during the evolution of dermatophytic and 
Onygenaceae species, providing flexibility for these species to 
occupy different ecological niches. 

Interestingly, LRT analysis in both M35 family and M36 family 
suggested that the duplicated Meps genes most likely underwent 
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Figure 6. Phylogenetic tree of 38 M36 family genes used for codon-based maximum likelihood analysis in PAML. Phylogenetic trees 
were collapsed with inconsistent nodes from different tree-building methods and poor statistical supports into polytomy. Branches a-j indicated 
putative duplication events in Onygenales fungi. The branches with significant evidence of positive selection are indicated as a thick branch. The 
putative positively selected residues along these branches were shaded in grey. 
doi:10.1371/journal.pone.0090225.g006 



positive selection during the evolution in dermatophytic and 
Coccidioides species. For dermatophytic fungi, significantly different 
selective pressures have acted on several lineages leading to both 
duplicated M35 family (branches a, e and J) (Figure 5) and M36 
family genes (branches a, b and g) (Figure 6). The results suggested 
that functional divergences of both M35 family and M36 family 
genes in dermatophytic fungi most likely occurred after gene 
duplication and before the speciation of these five dermatophytic 
species. The functional adaptations of these duplicated Meps in 
dermatophytic fungi might have helped these fungi survive in 
various niches because dermatophytes have many different 
environmental niches (soil, skin, hair), with different hosts (plants 
and animals), and with different cohabitating microbes. The 
characterization and identification of several Mep's functions in 
recent studies further indicated the important functions of Meps 
during host adaptation in dermatophytic fungi. For example, a 
43.5-kDa Meps belonging to the M35 family in the dermatophytic 



fungus M. canis was identified and it could digest keratin azure 
during infection of cats and guinea pigs [3]. However, the 
expression of genes encoding M36 fungalysins in A. benhamiae 
assayed by cDNA microarray analysis suggested that the M36 
family genes not only help A. benhamiae digesting keratin, but also 
facilitate the cutaneous infection of guinea pigs [54]. Both the M35 
family and M36 family genes have the capacity to degrade keratin, 
suggesting that the cooperation between the two family genes may 
have facilitated the fungi to adapt to various environmental niches. 

For Onygenaceae fungi, the positive selection pressure might 
have acted on the ancestral lineage (branch j) (Figure 5) leading to 
the three additional duplicated genes specifically in Coccidioides 
species, consistent with the results of our earlier study which 
suggested that significant selective pressure on the ancestral lineage 
of three additional duplicated M35 family genes from Coccidioides 
species may be associated with recent pathogenesis acquisition of 
M35 family genes in the Coccidioides species [13]. Noticeably, 
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although the M36 family genes from Coccidioides species and U. 
reesii were duplicated before the split of Arthrlodermataceae and 
Onygenaceae species, positive selection did not seem to have acted 
on the duplicated M36 family genes in Onygenaceae species. 
These results suggest that the M36 family genes in Coccidioides 
species and U. reesii are unlikely important functional genes during 
the evolution of these fungi. Of course, we can't exclude the 
possibility that these M36 family genes in Coccidioides and U. reesii 
may serve other yet unknown physiological functions. 

In our study, eleven and nine amino acid substitutions were 
predicted as positively selected sites respectively for the M35 and 
M36 family genes in Onygenales fungi (Table 2 and 3, Figure 5 
and 6). Most of the positively selected sites had substitutions of 
amino acids with different physical-chemically properties, such as 
differences in polarity, charge and/or hydrophobicity (Table. SI). 
It has been recognized that substitutions with chemically very 
different amino acids often have significant effects on the overall 
flexibility, structure stability and/or activity of proteins [55,56]. 
Therefore, although the explicit functional changes of the Meps 
brought by these positive selected residues identified in our study 
cannot be predicted at present, we propose that these residues may 
have allowed the development of new physiological functions of 
the duplicated Mep genes, providing valuable information for 
functional adaption of the duplicated Meps in dermatophytic fungi 
and Coccidioides species. Experimental investigations are needed in 
order test the functional effects of the potentially adaptive amino 
acid replacements discovered by our analysis. 

Conclusions 

In summary, our study analyzed the phylogenetic relationships 
among genes in both M35 and M36 families from Ascomycota 
fungi. Our results increase the current knowledge of the evolution 
of Meps in fungi and support the possibility that crucial 
physiological functions of both duplicated M35 family and M36 
family genes may have been developed during the evolution of 
fungi. Furthermore, our detection of significant positive selections 
in several lineages of dermatophytes and Coccidioides species 
indicates that rapid functional divergence of duplicated Meps 
have most likely occurred after gene duplications in these fungi. 
These species-specific adaptations could be important in host 
immune system interaction or in their survival in the environment. 
The hypotheses proposed in our analysis will provide valuable 
information for functional analysis of these genes in future studies. 
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